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In this work we introduce a low-energy Hamiltonian for single layer and bilayer black phosphorus 
that describes the electronic states at the vicinity of the gamma point. The model is based on a 
recently proposed tight-binding description for electron and hole bands close to the Fermi level. 
We calculate expressions for the Landau level spectrum as function of magnetic held and in the 
case of bilayer black phosphorus we investigate the effect of an external bias on the electronic band 
gap. The results showcase the highly anisotropic character of black phosphorus and in particular 
for bilayer BP, the presence of bias allows for a held-induced semiconductor-metal transition. 


In the last ten years the properties of crystals consist¬ 
ing of one or few atomic layers has been the focus of 
intense research. Such interest arose mainly due to the 
production of graphene in 2004, which has been shown 
to display remarkable electronic, optical and mechanical 
properties [l| . Since then, there has been a growing inter¬ 
est in the production of other low-dimensional crystals. 
The investigation of analogs of graphene has resulted in 
the discovery of several single layer crystals of different 
elements, such as Silicon (silicene) Q, Germanium (ger- 
manene) , as well as a class of materials known as tran¬ 
sition metal dichalcogenides Q . Some of these materials 
may soon find use in electronic devices, mainly due to 
the fact that in contrast with graphene, they present a 
band gap in their electronic spectrum, albeit with a lower 
carrier mobility. Among the most promising of these 2D 
materials is an allotrope of Phosphorus, known as black 
phosphorus (BP) [bMl l \ . which is that element’s most sta¬ 
ble crystal at room temperature and pressure. In bulk, 
BP is a narrow gap semiconductor with a orthorhombic 
structure that consists of atoms covalently bound into 
layers coupled by van der Waals interactions. Similarly 
to graphene, BP can be mechanically exfoliated to ob¬ 
tain samples with few or single layers, with the latter 
being known as phosphorene. The resulting material has 
a band gap that depends on the number of layers, vary¬ 
ing from 0.6 eV for five layers to 1.5 eV for a single layer, 
with carrier mobility in the range of ~ 1000 cm^ V“^s“^. 


The importance of a thorough understanding of the 
band structure and charge carrier dynamics in BP has 
led to a series of recent studies that obtained the elec¬ 
tronic dispersion using approaches such as first princi¬ 
ples calculations, k • p methods, as well as tight-binding 
models El- These calculations have shown evidence of 
a large anisotropy on the effective mass, as well as given 
estimates of the energy gap for single and multilayer BP. 
Calculations have shown the possibility of a topologial 
phase transition in few-layer BP, in which an external 
bias induces a band inversion [1^ . This would allow the 
development of devices in which the topological character 


of the material can be externally controlled. 

In this work, we consider the charge carrier dynamics 
in single layer and bilayer phosphorene by means of a 
continuum model obtained as the long wavelength limit 
of a recently proposed tight-binding model [l2| . In addi¬ 
tion to the anisotropy of the spectrum, another striking 
feature of the electronic bands obtained from this model 
is the hybrid nature of the electron and hole states close 
to the band edge in phosphorene, which display both a 
Schrodinger-like and Dirac-like character, which in turn 
is dependent on the direction of propagation. For the 
case of bilayer BP, we also consider the effect of an ex¬ 
ternal bias on the spectrum. We obtain results that show 
a bias-induced gap closure, which leads to the presence 
of zero-energy Landau levels. 

The paper is organized as follows: in section II we 
present the model Hamiltonian for single BP layers and 
analytical expressions for its Landau level spectrum. Sec¬ 
tion III extends that model for the case of the bilayer. Fi¬ 
nally, in section IV we present a discussion of the results 
and conclusions. 


SINGLE LAYER PHOSPHORENE 


The structure of each layer of BP has phosphorus 
atoms covalently coupled to three nearest neighbors. The 
resulting lattice resembles the honeycomb structure of 
graphene, however in phosphorene the sp^ hybridization 
of the 3s and 3p atomic orbitals creates ridges that result 
in a puckered surface (Fig . I). Using the tight-binding 
model proposed in Ref [l2| , we can write the Hamiltonian 
for single layer black phosphorus as 
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with eigenvectors given by [(I)a<Pb<Pd<Pc]^ where 
ua,b,c,d represent the on-site energies - which we hence- 
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FIG. 1: (Color online) Nearest neighbors in the phosphorene 
lattice. 


forth assume as equal to [/, with the A — C subscripts 
denoting the four sublattice labels shown in Fig. 1. The 
interaction terms are given in the appendix. By taking 
into account the symmetries of the phosphorene lattice, 
one can write a reduced two-band Hamiltonian for single 
layer black phosphorus at the vicinity of the Fermi level 
as 


nj _ ( U ^tAD^k) tAB{k) ^tAc{k)\ 

^ ~ \{tAB{k) + tAc{k)r U + tAD{k) ) 


which acts on the spinors 

il, = 1(^A + 4>d\ 


( 2 ) 

(3) 


From the Hamiltonian Eq.(2) one can obtain the energies 
for the bottom of the conduction band and the top of 
the valence band as Ec = 2ti + ^2 + 2ts + ts + 4^4, and 
Ey = — (2ti + t 2 + 2ts + ts) + 4^4. That leads to a gap of 
A 1.52 eV. 

By diagonalizing the Hamiltonian Eq. (2) one can ob¬ 
tain the following dispersions: 


E{kx^ky) = U E At 4^0.0^ {kxdi) co^{kyd2) 


±|4[t^ + tg + 2tits cos { 2 kyd 2 )] cos^ (kxdi) 

+ [^2 + ^5 + 2 t 2 t^ cos { 2 kyd 2 y\ 

+4^3 [^2 cos {kyd 2 ) + ts cos {3kyd2)] cos (kxdi) 

. 1/2 

+4ti[t2 + h] cos (kxdi) cos {kyd 2 ) > , (4) 

where di = aismai/2 and d 2 = ai cos(ai/2 + a 2 cos/3, 
with the positive (negative) sign corresponding to the 
conductance (valence) band. Eigure 2 shows a plot of Eq. 
(4) centered at the gamma point (black lines), where the 
strong anisotropy of the spectrum is evident. 

A simple calculation shows that the eigenstates of the 
Hamiltonian Eq.(l) can be found as 

= M > (5) 

V2 j 

where A = ±1, with the same sign convention as Eq. (4) 
and 

Ok = tan“^(C/T^) (6) 


where 

C = —2ti cos (/ca,di) sin (/c^ai cos(ai/2) + t 2 sin (/c^a 2 cos/3) 
+ 2^3 cos {kxdi) sin [ky{ai cos((ai/2) + 2a2 cos/3)] 

—ts sin [ky{2ai cos(ai/2) + a 2 cos/3)] (7) 


and 

D = 2ticos{kxdi) cos{kyaicosai/2) E t2Cos{kya2Cos(3) 
E2ts cos (kxdi) cos [ky{ai cos(ai/2) + 2a2 cos/3)] 

Eh cos [ky (2ai cos((ai /2) E a 2 cos /3)]. (8) 

Although these results show some similarity to the re¬ 
sults for graphene it can be seen that for phosphorene 
the phase angle does not correspond to the polar angle 
of the momentum vector. 

Continuum approximation: By expanding the struc¬ 
ture factors around k = 0 (Gamma point) and retaining 
the terms up to second-order in /c, one can write a long- 
wavelenght approximation for the Hamiltonian Eq.(2) as 

H = ( ^ 

^ \6 + jAl + 'lykl - ixky Uo + ijAl + Vyk^ J 

(9) 

where 


r]x = -2t4[aisin((ai/2)]^, 

rjy = — 2 t 4 [ai cos((ai/2) + a 2 cos/3]^, 

lx = -(G+t3)[aisin((ai/2)]^ 

7 y = —ti[ai cos((ai/2)]^ — t 3 [ai cos((ai/2) + 2 a 2 cos/3]^ 
—t 2 {a 2 coSyd)^/2 — h[2ai cos((ai/2) + 02 coSyd]^/2, 
d = t2 E t^ E 2(ti E ts)^ 

Uo = 4^4, 
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FIG. 2: (Color online) Low-energy dispersion of phosphorene 
from tight-binding model (solid black lines) and continuum 
approximation (dashed red lines). 


X = t 2 a 2 cos/3 + 2 t 3 [ai cos(ai/2) + 2 a 2 cos/3] — 

t^{2ai cos[ai/2) + a 2 cos/3] — 2tiai cos(ai/2). (10) 

By substituting the hopping parameters in the above 
expressions we obtain the following values: uq = —0.42 
eV, rja, = 0.58 eV- Vy = 1-01 eV- A^, (5 = 0.76 eV, 
X = 5.25 eV- A, 7 ^ = 3.93 eV- AA and jy = 3.83 eV- AA 
The eigenvectors are [4>i with the <pi ^2 spinor 

components now corresponding to envelope functions as¬ 
sociated with linear combinations of the amplitudes for 
each sublattice site. The form of Hamiltonian Eq. (9) 
is similar to the one presented in Ref. [7], which was 
obtained within the k • p approximation with parame¬ 
ters chosen in order to fit the band structure obtained 
from first principle calculations. In the present case, 
however, the parameters include the contribution from 
different hopping terms between neighboring lattice sites, 
as well as the lattice geometry, and thus can be under¬ 
stood as presenting a direct link between the microscopic 
tight-binding description and the continuum approxima¬ 
tion. Moreover, within this model additional momentum- 
dependent terms can be added to better approximate the 
spectrum at higher energies by including higher-order k 
terms in the structure factor expansion. Dispersion rela¬ 
tions for electrons and holes are then given by 

E = uo + Vxkl + Vykl ± 

( 11 ) 

where the plus (minus) sign yields the conduction (va¬ 
lence) band. As shown in Fig. 2, there is good agreement 
between the continuum and the tight-binding results for 
energies in the range —2.0 to 1.5 eV. It can be seen, 
from the spectrum of Eq. (11) that, although BP has 
an anisotropic dispersion, it does not correspond exactly 
to the spectrum of a simple anisotropic system with an 
elliptical dispersion, due to the additional term propor¬ 
tional to ill radical. However, as shown below. 


for low energies a simple anisotropy on the effective mass 
can be recovered as an approximation. 

Effective masses: From the spectrum Eq. (11) one can 
estimate the effective masses of electrons and holes in 
BP. Taking into account the anisotropy of the system, 
one can readily find, for the x direction: 


e h 

2 ( 7 /^+ 7 ^)’ ^ 2 {xx-Vx)' 


( 12 ) 


For rriy one finds, for small values of ky^ 


2(7?y ± 7 j^ ± xV2(5) ’ 


(13) 


where the plus (minus) sign corresponds to electrons 
(holes). The resulting effective masses are m% = 0.846 
mo and m^ = 1.14 mo, m^ = 0.166 mo and rriy = 0.182 
mo, with mo being the mass of a free electron. In com¬ 
parison, the values of effective masses quoted in Ref. [8] 
are m| = 0.7 mo and m^ = 1.0 mo, and rriy = rriy = 0.15 
mo (in that case, the choices of x and y labels were the 
opposite of ours). One then can use these masses to ob¬ 
tain an approximation for the spectrum Eq. (11) as (for 
electrons): 


E={uq^5)^ 


2m% 


—e 

2ml y' 


(14) 


and a corresponding expression for holes. 

Eigenstates: The continuum approximation Hamilto¬ 
nian Eq.(9) can be rewritten in a more compact form 
as 


where 


H 


ei e2e*®''\ 
€ 26 “®^'= €1 ) 


(15) 


= = + (16) 


Ok = tan ^[2xky/{f+ - /_)], (17) 


where we defined 


f± = (wo ± (5) + {rjx ± lx)kl + {rjy ± Xy)kl, (18) 

where, for = 0, the /+ and /_ expressions yield the 
dispersions for the conduction and valence bands, respec¬ 
tively. Thus, using this polar notation, one can readily 
obtain the eigenstates as 

^a = T| I , (19) 

V2 j 

where A = ±1, with the positive (negative) signs cor¬ 
respond to electrons (holes). These expressions are for¬ 
mally similar to the states of Eq. (5), which are valid 












for the whole Brillouin Zone and, as before, the angle 
Ok does not correspond necessarily to the polar angle as¬ 
sociated with the momentum vector. In fact, since the 
denominator in Eq. (18) depends only on even powers of 
the momentum components, the polar angle will assume 
values in the range —Oc < Ok < Oc^ where < 7r/2 is an 
energy-dependent critical value corresponding to kx = 0. 
From the form of Eq. (18) it is seen that as the energy in¬ 
creases Oc approaches a maximum value and then decays 
to zero. One consequence of that behavior is the fact 
that, although a pseudospin may be defined for charge 
carriers in phosphorene for a certain energy range, the 
Berry phase is nevert her less zero, due to the vanishing of 
the winding number around the E point. 

Landau levels: In order to calculate the Landau lev¬ 
els for phosphorene, let us consider the Hamiltonian Eq. 
(9) with a magnetic field perpendicular to the plane of 
the layer, and use the gauge A = {—By,0,0) and the 
substitution k —iV. Since the Hamiltonian does not 
depend on x, we can assume 0 i, 2 (^,^) = 
with 01 = (0^ -h 0 d)/ 2 and 02 = (0 b + 0c)/2. Thus we 
obtain the following pair of coupled differential equations 

d? 

K + Vx{kx + (iyf - 

+ [<5 + lx{kx + Pyf - = E(j)i 

(f 

[uo + r]x{kx + Pyf - 

+ [<5 + ^x{kx + ~ ~ E(j)2(20) 

where 0 = eB/h = , with £b being the magnetic 

length. Let us now set kx = 0 without loss of generality 
and rewrite the Hamiltonian in terms of ladder operators, 
acting on the spinor components 0± = (0i ± 02 )/\/ 2 , 



Thus, we can readily obtain a Hamiltonian in terms of 
the a operators as 


( 22 ) 

where 1 is the unit matrix, ax and az are Pauli matrices 
and 



FIG. 3: (Color online) Landau levels as function of magnetic 
field. The linear approximation is shown as the red dashed 
curves. 

and Ay = {jx ~ ly)I A plot of the Landau levels as 
function of magnetic field is shown (black dots) in Fig. 
3, for electrons. 

Although the actual spectrum deviates from the linear 
dependence on magnetic field for large fields, for B < 
30 T the spectrum still shows an approximately linear 
dependence. In this regime, one can obtain an expression 
for the Landau levels by means of the following ansatz: 

This ansatz can be justified by the fact that its introduc¬ 
tion leads to an approximate Hamiltonian in which an 
additional term proportional to y is added to the y-mass 
term (see, e.g. Eq. (13)). Thus, using the above ansatz 
allows us to obtain a block diagonal Hamiltonian where 
the block corresponding to the electron branches is 

= Uq + (5 + 277 + 0 ( 0 ^^ 0 ^ -j- 1/2) -|- A+0(of^Of^ + OfOf) 

y 2 

——I3{a^a^ + aa — aa^—a^a). (26) 

We now define 

/'i = = A+ - (27) 

that allows us to rewrite the Hamiltonian Eq. (26) as 
He = uo-\-6 2jaif3{a^a + 1 / 2 ) - 1 - /i 20 (o^^o^^ + o^o^)- (28) 
Next, one can perform a Bogoliubov transformation 


— uq “h ^ T 277+0 (o^'^o^ T 1/2) A+ 0 (o^'^o^'^ -|- 0^0^), (23) 


c = wa -h va 




ci = wa^ 


(29) 


and 

£- = uq — S -\-2r]_P{a^a1/2) A_(3{a^-\-aa), (24) 

where we defined 77 + = 77 ± 7 and A+ = A 77 ± Ay, with 
0 ~ {Vx + Vy)/‘^i 7 ~ {'Jx + ly)l‘^^ ~ iVx ~ Vy)/‘^ 


with = 1, for which to w = coshz^, v = sinhz^, 

tanh2z/ = 1 ^ 2 /yi. That gives us 


w = 


Li 


\f2 L - nl 


+1 


1/2 


V = 


Li 


nl/2 


\/2 1 -++ 


-1 


( 30 ) 
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Finally, one can readily obtain the transformed Hamilto¬ 
nian for the electronic branches as 

l~Le =(^ + 1^0 + C -j- 1/2), (31) 

where 

lOq = (32) 

and (3 = eB /A similar approach yields, for the 
hole block, 

Bh = —S -i- uq — hujh{d^d + 1/2), (33) 

where the d operators are obtained from the a ladder 
operators via another Bogoliubov transformation, and 
where 

~ ^2 5 (34) 



The spectra obtained from Eq. (31) is shown as dashed 
red lines in Fig. 3 for Landau indices n = 0, to 6 . Similar 
expressions for the Landau levels in single layer phospho- 
rene where obtained recently by means of a perturbative 
calculation in ref. [l3|, which was based on the same 
tight-binding model employed here. However, in contrast 
with these results the present approach can be readily 
generalized for the bilayer case, as we show below. 


BILAYER PHOSPHORENE 


For the case of two coupled phosphorene layers, one 
now has to consider 8 sublattices, which we label 
A, H,C, D for the lower layer and and D' for 

the upper one. Using the tight-binding model of ref. [l2| 
one obtains the following Hamiltonian 


Bk 


[Hi H,\ 


(36) 


acting on the spinors T = 

[(j>A 4 >b <t>D (f>c (P'a 4>'b ^'d where 

the Hi ^2 blocks contain the interaction terms connecting 
sublattice sites within the same layer. 


/ Ul ,2 tAB{k) tAD{k) tAc{k)\ 
tAB{ky rti,2 tAc{ky tAD{k) 

tAD{k) tAc{k) Ui^2 tAB{k) 

\tAc{ky tAD{k) tAB{ky Ui ^2 ) 


(37) 


sites located in adjacent layers; here, these correspond 
to the sublattice sites A, H, C' and D' with [He] 13 = 
tAD'{k), [Bc]i 4 = tAC'{k), [Bc ]23 = tBD'{k) = tAC'{ky 
and [He] 24 = tBC'{k) = tAD'{k), with the remaining ele¬ 
ments being zero. The expressions for the coupling terms 
are given in the appendix. In the continuum approxima¬ 
tion, the coupling terms become 


tAB{k) = 5i+'jikl + 'j2kl+ixiky, 
tAc{k) = 52+x?.kl+iX2ky, 
tAD{k) = 53 + riikl+r] 2 ky, 
tAD'{k) = 5i+ rjskl + tjAI, 
tAC’{k) = 55 + XAkl + xAl + ixzky (38) 


where (5i = -2.85 eV, <52 = 3.61 eV, 83 = -0.42 eV, 
^4 =-0.06 eV, 85 = 0.41 eV, rji = 0.58 eV- A^, r /2 = 1.01 
eV- A^, 7 i = 3.91 eV- A^, 72 = 4.41 eV- A^, 73 = —0.58 
eV- A^, xi = 2.41 eV- A,X 2 = 2.84 eV- A, r )3 = 3.31 eV- 
A^, r ]4 = 0.14 eV- A^, 74 = —0.56 eV- A^, 75 = 1.08 eV- 
A^, and X 3 = 1.09 eV- A. 

The above Hamiltonian leads to a system of 8 cou¬ 
pled equations. However, as we show below, one can still 
obtain approximate analytical solutions. The eigenvalue 


(39) 


equation can be rewritten as 



(Hp H'- 

) = E-^' 

where 



II 

+ H2- \H3 

K ^#1 

—i-y 1 

H 0 +H 2 + \H3 

Hm = 

[Ho -H 2 - \H3 
1 

— 

H 0 -H 2 + \H3 


(40) 


(41) 


and 


h: = 


-hH3 


0 


(42) 


where 1 is the 2 x 2 unit matrix, A denotes ui — U 2 and 
we assume U 2 = —ui^ and 


Ho = 


i ° 

\tAB{k)* 


tAB{k)\ 

0 


(43) 


and 


TT _ (tAD{k) tAc{k)\ 

^ \tAc{ky tAoik) J ' 


f tAD'{k) tAC'{k)\ 

\tAC'{ky tAD'{k)) 


(44) 


(45) 


and the eigenvectors are the 8 -component spinor T' = 
[^pp ^rnp ^pm which the four sets of 2 - 

component spinors are 


with ui ^2 being the onsite energies for upper ( 1 ) and lower 
(2) layers. The He blocks contain the couplings between 


/ _ 1 f(f^A ^ ^ 4^A' + 4d'\ 

2 a /2 \4b + 0 c + 4b' + 4c') 












FIG. 4: (Color online) Band structure of bilayer black phos¬ 
phorus at the vicinity of the F point, obtained form a tight- 
binding model (black solid lines) and the continuum approach 
(blue circles). 


I _ 1 ((j^A — (j^D — (j^A' + 

- 0c - + 0C' ) 


I _ ^ (4>a + 4>d- - 4>d'\ 

2>/2 (t^c — 0 S' - 0 C' / 




^ f<pA — <pD + 4 >A' — 0 dA 

2\/2 \4>b - 0 c + 0 S' - 0 C ') ’ 


(46) 


A further approximation can be made by taking into ac¬ 
count the fact that the off-diagonal blocks give rise to 
a small perturbation to the spectrum and can thus be ne¬ 
glected in a first approximation, leading to the following 
pair of eigenvalue equations: 


(Ho + H2 + \h^-e -zfi WV’ppA n 

V i|l Ho + H2-\H:,-E)yPmp) ’ 

(47) 

and 

Fo - ^2 + -ifl 

i|l Ho-H 2-^H3-EJ y^Pram) 

(48) 

In this case, by solving Eq. (48) one obtains the 4 in¬ 
ner families of branches (i.e. closer to the Fermi energy) 
whereas Eq. (49) leads to the outer families of levels. 

In the absence of bias, these systems of equations are 
reduced to 4 copies of Eq. (9) although with different 
parameters. The resulting 8 bands are labeled 
and the parameters corresponding to the four low-energy 
branches are shown in Table I, with the effective masses 
given as multiples of the electron mass mo, with indices 
in decreasing order of energy. Eor finite bias, the Hamil¬ 
tonians Eq. (47) and (48) can be diagonalized. In order 


FIG. 5: (Golor online) Low energy espectrum of bilayer phos- 
phorene as function of the energy difference between layers 
obt; 

Har 



FIG. 6: Landau levels as function of magnetic field, for Ui = 
0.74 eV and U 2 = -Ui. 


to show that, let us first recall Eq. (16) and rewrite the 
2x2 diagonal blocks in Eq. (47) as 


H0 + H2 + 

and 

Ho+H2-^H3 



e" e' 2 'e<\ 
e' 2 'e-*®^' e" ) 


(49) 


(50) 


with the e'l 2 , 2 polar angles are defined as in 

Eqs. (17)-(19). Thus, after some straightforward alge¬ 
bra, we can obtain the four energy bands associated with 
Eq. (47) as the solutions of the equation 

[{E - - {e',n{E - e'lf - (e")"] = -(|)" 
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-0,38 


> 

W 



-0,42 k 


-0,44 k 


with 


/i = (58) 

^2 

and the other terms defined as before. It can be easily 
seen that, as A ^ 0 we obtain a ±1, 6, c ^ 0, as 
expected. For the valence band, the result is similar, 
with 



(59) 


where 


FIG. 7: Landau levels as function of on-site energy, for B = 10 
T, and U 2 = —Ui. 




■ i—b 


(60) 





COS 


{0', - 01) + {E- e'^){E 


e'/) . (51) 


For the range of energy and momenta of interest, one can 
safely assume cos — O'l) 1. In that case, Eq. (51) 
becomes 

[{E-e\-e'^){E-e'i-4)-[^)^] 

x[{E - e'l + e')(E - e'i + e") - (|^)'] = 0. (52) 

One can then obtain expressions for the energies of the 
low-energy bands at the F point as function of A as 

= + +(y) 

E, = -(5i - ^2 + ^3 + +{^- (53) 

Eigenstates: Plane-wave eigenstates for the inner 
bands can be obtained from the Hamiltonian Eq. (47) 
as, for the conduction band: 


^c{k) = A, 



(54) 


with 


dc 


{E - e'l) 



(55) 


and 


Cy 


{E-e'O, 



2 [{E - e'l)^ - ef + A‘^h'/4] 
A [E-e'{ + h'{E-e[)] 


(61) 


(62) 


where h' = 1/h. The normalizing constants are given by 

^c,v = [1 + \dc,v\^ + • 

Eigure 4 shows the spectrum of bilayer BP obtained 
from the tight-binding model (black solid lines) and con¬ 
tinuum approaches (blue circles). As in the case of the 
single layer, the continuum results show a good agree¬ 
ment with the tight-binding data for the range —1.5 to 
1.5 eV. The effect of biasing on the gap is shown in Eig. 
5 with data obtained from both the original 8x8 tight- 
binding Hamiltonian (black solid lines) as well as from 
the analytical expression Eq.(54) (blue dashed lines). 
The results show a good agreement, with a discrepancy 
of ~ 4 meV. Eor values of A above ^1.5 eV, the conduc¬ 
tion and valence bands overlap, and the system becomes 
metallic. 

Landau levels The equations above lead to a set of 4 
electron and 4 hole families of Landau level branches. In 
the absence of biasing (i.e. A = 0), Eqs. (48) and (49) 
can be solved analytically in a similar fashion as in the 
case of single layer, with the parameters modified by the 
presence of interlayer coupling. Thus, the expressions for 
the different families of Landau level branches have the 
form 


and 


2 [{E - e[)^ - e'i ^ A^h/4] 
A [E- e[ + h{E - e'/)] 


(56) 


( 57 ) 


E = 5e±hu}e{n+ 1/2), (63) 


where uji = the i indices denote different 

combinations of the coupling terms, with the positive 
sign corresponding to frequencies of electron branches 
{£ = and the negative sign is assigned to the 
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TABLE I: 


£ 

Hi 

iv 

V 

vi 

Se 

0.515 eV 

0.165 eV 

-0.947 eV 

-1.413 eV 

mi 

0.65 

1.23 

0.72 

2.74 

IJly 

0.17 

0.16 

0.17 

0.18 


hole branches {i = The values of 5i and the 

effective masses are displayed in Table 1. 

In the presence of an external bias, a numerical ap¬ 
proach becomes necessary. Figure 6 shows the depen¬ 
dence of the energy levels on the magnetic field, for a 
finite bias (A = 1.48 eV). In this case, for the range 2 T 
^ B ^ 12 T the n = 0 LL becomes doubly degenerate 
and weakly dependent on 5, a situation that is analo¬ 
gous to the case of single layer graphene. This analogy is 
reinforced by the fact that the remaining levels become 
unevenly spaced. 

The dependence of the Landau levels on the bias, for 
a fixed magnetic field, is shown in Fig. 7. It is seen that 
the presence of the external electric field tends to close 
the gap for a certain critical value of the bias. Moreover, 
the branches tend to become degenerate. 

CONCLUSIONS 

We have presented a continuum description of single 
layer and bilayer black phosphorus, starting from a tight 
binding model that reproduces the results of first prin¬ 
ciples calculations. Using this model we obtained the 
spectra of electrons and holes at the vicinity of the Fermi 
level at the gamma point and calculated the Landau level 
spectrum for both systems. For the case of bilayer BP 
we considered the effect of interlayer bias by introduc¬ 
ing a layer-dependent on-site energy in the model. This 
showed that the presence of bias can close the electronic 
band gap, converting the material into a metal for a crit¬ 
ical value of on-site energy difference. Correspondingly, 
the Landau level spectrum shows the appearance of dou¬ 
bly degenerate branches with a zero energy level weakly 
dependent of magnetic field for on-site energies above the 
critical value. This result agrees with recent ab initio cal¬ 
culations for few-layers black phosphorus 0 and can be 
exploited as the basis for future gate-tunable electronic 
devices. Furthermore, by taking into account additional 
interlayer hopping terms in a tight-binding description, 
the present model can be readily extended to deal with 
multilayer black phosphorus. 
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APPENDIX 

The structure factors obtained from the tight-binding 
model of ref. 0 are given by the expressions 

tAB{k) = 2ti COS (/ca,ai sin((ai/2))x 
ex.p[—ikyai cos {ai/2)] 

-\-2ts cos {kxai sin((ai/2)) x 
ex.p[iky{ai cos (<ai/2) + 2a2 cos /3)] (64) 

tAc{k) = t 2 exp[i/c^a 2 cos/3] 

+^5 ex.-p[—iky{2ai cos((ai/2) + 02 cos/3[65) 

tAD{k) = 4^4 cos (/ca,ai sin((ai/2)) 

X cos[ky{ai cos((ai/2) + 02 cos(/3)], (66) 

where ai is the distance between nearest neighbor sites 
in sublattices A and B or C and D, and 02 is the distance 
for n.n. sites of A and C or B and D; ti and t 2 are the 
corresponding hopping parameters for nearest-neighbor 
couplings. Due to the symmetry of the lattice, we have 
that t^jjik) = {t'j^^{k)Y^ t'cB{k) = tAoik)^ ^'BD^k) = 
{t'Ac{k)Y^ and tBc{k) = tAoik). The bond angles are 
shown in Fig. 1, and the parameters are ai = 96, 5°, 0^2 = 
101,9°, cos/3 = — cos 0 ^ 2 /cos (a 1 ai = 2.22 A, 02 = 2.24 
A. The hopping parameters are, in eV, ti = —1.220, t 2 = 
3.665, ts = -0.205, U - 0.105, and h. - 0.055 

For the case of bilayer BP, the additional coupling 
terms are 

tAD'{k) = 4t^ COS (/ca,2ai sin((ai/2)) X 

cos {ky{ai sin((ai/2) + 02 cos/3)) 

-\- 2 t 2 cos {ky{ai sin((ai/2) + 02 cos/3),(67) 

and 

tAC'{k) = 2t^ cos (/ca,ai sin((ai/2)) exp[i/c^a 2 cosyd] 
cos {kxai sin((ai/2)) 

X ex.-p[—iky{2ai sin((ai/2) + 02 cos/3)],(68) 

with ti = 0.295 eV, = 0.273 eV, = —0.151 eV and 
ti = -0.091 eV. 
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